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Abstract 

We analyze the FEAST method for computing selected eigenvalues and eigenvectors of large sparse matrix 
pencils. After establishing the close connection between FEAST and the well-known Rayleigh-Ritz method, 
we identify several critical issues that influence convergence and accuracy of the solver: the choice of the 
starting vector space, the stopping criterion, how the inner linear systems impact the quality of the solution, 
and the use of FEAST for computing eigenpairs from multiple intervals. We complement the study with 
numerical examples, and hint at possible improvements to overcome the existing problems. 
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1. Introduction 

In 2009, Polizzi introduced the FEAST solver for generalized Hermitian definite eigenproblems PQ. FEAST 
was conceived as an algorithm for electronic structure calculations, and then evolved into a general purpose 
solver. In this paper we describe the mathematical structure of the algorithm and conduct an investigation 
of its robustness and accuracy. 

FEAST belongs to a family of iterative solvers based on the contour integration of a density-matrix 
representation of quantum mechanics; the result of the integration is a subspace projector that plays a central 
role in a Rayleigh-Ritz method. The prominent member in this family was developed by Sakurai and Sugiura 
in 2003 |2j; this work, which inspired a number of generalizations and high-performance implementations [3J 
H], targets non-Hermitian eigenproblems. By contrast, FEAST promises to deliver performance on sparse 
Hermitian problems. Such problems can be solved by a number of alternative packages like ARPACK [S] 
and TRLan [5] and the solver implemented in PARSEC [HIH]. 

Since in ab-initio electronic calculations typically one is interested in the lowest part of the eigenspectrum, 
we investigate FEAST 's behavior for the computation of a subset of eigenpairs lying inside a given interval. 
Additionally, we study its strengths and weaknesses when a large portion or all of the spectrum is sought 
after. In our analysis, we use a number of matrices from practical applications. From our experiments, we 
found that while for specific scenarios FEAST is accurate and reliable, in general it lacks robustness. 

Our analysis touches upon three main features of the solver: 1) critical input parameters, 2) the stopping 
criterion, and 3) the quaUty of the results. 
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1) In addition to a search interval, FEAST's interface requires the user to specify the number of eigenval- 
ues present within the interval. Although in some ab-initio simulations this number can be accurately 
estimated, in general it is not possible to obtain it cheaply; since the completeness of the computed 
eigenpairs greatly depends on it, this initial guess is critical. The user can also specify a starting vector 
base which FEAST uses to initialize the solver: while by default FEAST uses a random set of vectors, 
the convergence rate of the solver greatly depends on the actual choice. We elaborate on how different 
starting bases affect convergence speed and robustness. 

2) The original implementation of FEAST employs a stopping criterion based on monitoring relative 
changes in the sum of all computed Ritz values. We identify cases where this criterion does not reflect 
the actual convergence and propose an alternative criterion based on per-eigenpair residuals. 

3) We analyze the quality of the solution computed by FEAST by means of residuals and orthogonality. 
On the one hand we found that the achievable residuals are affected by the accuracy of the linear solver 
used within the algorithm. On the other hand, we distinguish between local and global orthogonality, 
i.e., orthogonality among eigenvectors that were computed with a single interval, and among vectors 
from separate intervals, respectively. While the local orthogonality is guaranteed by the Rayleigh-Ritz 
method, depending on the spectrum of the eigenproblem the global orthogonality might suffer. 

The paper is organized as follows. In Section |2] we illustrate the underlying mathematical structure of 
the algorithm. Section [3] contains experiments and analysis concerning the different aspects of the solver. 
Section ^ examines the suitability of FEAST as a building-block for a general purpose solver, working 
on multiple intervals to compute a larger portion of the spectrum. We conclude in Section [5] suggesting 
improvements to broaden FEAST's applicability. 

2. FEAST and the Rayleigh-Ritz method 

Let us consider the generalized eigenproblem Ax = XBx, with Hermitian matrix A G C"^" and Hermitian 
positive definite B S C"^". The objective is to compute the eigenpairs whose eigenvalues lie in a given 
interval I\ ~ [A, A]. Since FEAST is an instantiation of the Rayleigh-Ritz method, we start with a short 
review. 

In the following, in order to denote an eigenpair (A, x) with A G /a, we will sloppily say that the eigenpair 
is in the interval Ix- 

2.1. Rayleigh-Ritz theorem 

We present the Rayleigh-Ritz method in its orthonormal version, which ensures a minimal residual. The 
method relies on the following theorem. 

Theorem 2.1 (Rayleigh Ritz, [9]) 

Let U be a subspace containing an eigenspace X <Z lA of the generalized eigenproblem Ax — XBx. Let U 
be a basis of vectors for U = span([/), a left inverse of U, and Aij — AU, Bij — BU the so-called 
Rayleigh quotients for A and B. If (A, W) (A = diag(Ai, X2, ■ ■ ■)) are primitive Ritz pairs of the reduced 
problem, i. e., AuW = BuWA, then (A, X) are Ritz pairs for the original eigenproblem with X = UW and 
span(X) = X . 

In a neighborhood of the exact solutions for the primitive Ritz pair, we expect that the Rayleigh-Ritz 
theorem also applies approximately, leading to the following strategy. 

1. Find a suitable basis U ior U. 

2. Compute the Rayleigh quotients Au = U^AU, Bu = U^BU. 

3. Compute the primitive Ritz pairs (A, W) of AuW = BuWA. 

4. Return the approximate Ritz pairs (A, UW) of AX = BXA. 
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5. Check convergence criterion; if not satisfied, go back to Step 1. 

Let us point out that obtaining an accurate approximation of (A,X) is not an obvious consequence of 
computing primitive Ritz pairs. One must ensure that both Ritz values and the corresponding Ritz vectors 
converge to the desired eigenpairs. Conversely, it must hold that each eigenpair in I\ corresponds to a 
primitive Ritz pair in the same interval; see 



2.2. The algorithm 

The FEAST algorithm implements the above Rayleigh-Ritz method, with a particular choice for computing 
U: 

dz{zB - Ay^BY, (1) 

2m Jc 

where C is a curve in the complex plane enclosing the selected interval I\. The expression {zB — A)^^B 
is normally referred to as the eigenproblem's resolvent; Formula ([T]) can be interpreted as the projection of 
the set of vectors Y onto a subspace U containing the eigenspace. Pseudo-code for FEAST is provided in 
Algorithm [T] 

Algorithm 1 Skeleton of the FEAST algorithm 

Input: An interval Ix = [A, A] and an estimate M of the number of eigenvalues in I\. 
Output: M < M eigenpairs in Ix- 

1: Choose Y e C"^*^ of rank M and compute J7 := 2?? /c (^^ " A)-^B Y; 

2: Form the Rayleigh quotients Au := U^AU, Buj= U^BU^ 

3: Solve the size-M generalized eigenproblem AuW = BijWK; 

4: Compute the approximate Ritz pairs (A,X :— U ■ W); 

5: If convergence is not reached then go to Step[l) with Y := X. 



Having established the close connection between FEAST and the Rayleigh-Ritz method, the next section 
is devoted to the theoretical framework for the resolvent and the contour integration ([!]) . 

2.3. Integrating the resolvent 

In this section we define the concept of resolvent and illustrate its functionality within the Rayleigh-Ritz 
method. The objective is to show that the subspace U is approximated by the integral operator in Equa- 
tion Q. We recall that a generalized Hermitian definite eigenproblem has n real eigenvalues Ai, . . . , A„ 
and i3-orthonormal eigenvectors Xi, . . . , x^- 

Let us consider a single eigenpair (A^, Xk), and let z G C. For z ^ £ {!,..., n}, 

B^^{zB - A)xk ^ (z - Xk)xk, 

and 

{zB - A)-^Bxk = {B-\zB - A))-^Xk = {z - X^T^Xk. (2) 
Define the resolvent operator G{z) as 

G{z) {zB~A)-^B 

and let be a closed curve in the complex plane enclosing only the eigenvalue A^,. Thus, the integral 
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dz G{z)xk 

equals the residue of G{z)xk localized at the pole in Afe. Using Equation ([2|, we obtain 

1 /" \ f dz 1 

— / dz {zB - Ay^ Bxk ^ / —Xk^^r—2'n:ixk = Xk. (3) 

2m Jr.. 2m Jr^ z ~ \k 2m 



By contrast, the residue around any other pole returns 0. 

Now let C be a curve enclosing a subset {A^ : k £ 1} of the eigenvalues. Combining Equations ([s]) for 
k G I and splitting the path integral over C into a sum of integrals over closed curves Ck containing just one 
eigenvalue A^. each, one obtains 

^ / d. = E i / = E 4.x, = { J- ^^^^^ . (4) 

So far we have shown how the projection operator acts on the full space of eigenvectors. This is a well 
known property of the resolvent of an eigenproblem |10l Ch. 3]. This property can also be visualized by 
combining the B-orthogonal eigenvectors Xk, k £ I, into an rt-by-|/| matrix X, and by comparing Q with 
the application of the projector 



Q = XX" B = Y^ Xkxf B with Q2 = g 



fee/ 



to an eigenvector Xj, 



^x, = 2^{xux'^)Bx, = E - 1 o: otherwise 

fee/ fee/ ^ 

One then concludes that the operators 2^ dz G{z) and Q, when applied to a set of eigenvectors, produce 
the same results. 

Let Y = {yi,y2, Vm}, then 



a J el 



— [ dz izB - A)-^ BY =QY = XX^BY 

271-1 Jr 



projects each yj onto the eigenspace X — span(X)j^ In this sense, the matrix U computed in Algorithm [T] is 
a reasonable attempt to fulfill the requirements of the Rayleigh-Ritz theorem. In practice, the integral (llj) 
must be evaluated numerically, using a scheme such as Gaufi-Legendre; for details we refer to the original 
publication [T1 and to the literature on numerical integration, e.g., 
A numerical integration scheme leads to an approximation 



- Vwfe(zfeS- BY 



U - „ 
2m 

fe=i 

where the points Zk lie on the curve C. The accuracy of the approximation, as well as the computational 
complexity, are determined by the number of integration points. In practice, Gaussian rules with 7-10 nodes 
already achieve satisfactory results. For each integration point, a linear system {zkB — A)Uk = BY of size 
N with M right-hand sides must be solved. 

If the curve C is chosen to be symmetric with respect to the real axis (e.g., a circle or ellipse), then 
considerable computational savings are possible. In this case the numerical integration must cover only the 
half-curve in the upper complex half-plane due to the symmetry G{z) — G"{z) [IJ. 

3. Analysis and experiments for the FEAST algorithm 

In this section we discuss the issues arising when employing FEAST to compute the eigenpairs in a single 
interval: size and choice of the search space, stopping criteria and impact of the linear solver on the accuracy 
of the solution. 



* In Polizzi's original paper, U equals XX^Y instead of XX^ BY . This is correct only when Y is chosen to be a random 
set of vectors, since BY does not alter the random nature of Y. Thus, this equality is only valid in the first iteration of the 
solver. 
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3.1. Size of the search space U 

Section [23] shows that the size of the eigenspace X is determined by the number M of columns of U which in 
turn corresponds to the number of columns of Y . This number is equivalent to the number M of eigenvalues 
(counting multiplicities) lying in the given interval I\. In practice, M is not known a priori, and an estimate 
M has to be used instead. In the following, we discuss the consequences of the cases M > M and M < M. 

Case M > M. If the estimate is largerj±ian the actual number of eigenvalues in I\, then the spectral 
projector Q has rank M, and thus the n-hy-M matrix U = QY is rank-deficient and does not have a full rank 
left inverse. As a direct consequence, Bjj is rank deficient, and the eigenpairs of the reduced eigenproblem 
might be incomplete and not necessarily i?-orthogonal. 

In the original implementation of FEAST, the rank deficiency is detected by measuring the "positive 
definiteness" of Bu through a Cholesky decomposition. A drawback of this approach is that it does not 
indicate which of the columns of U are linearly independent, forcing one to select a subset arbitrarily. This, 
in turn, might result in the presence of spurious eigenpairs. A reliable, but more expensive, approach consists 
of computing an SVD or a (rank-revealing) QR decomposition of the matrix U [E]. Since the rank of U 
equals the number of eigenvalues in Ix, such a rank-revealing decomposition can be used to safely restart 
the process with Y = UX, where U includes the linearly independent columns of U. 

Case M < M. The space spanned by U does not contain the whole eigenspace corresponding to the 
eigenvalues within the integration contour, and therefore the vectors generated by integrating the resolvent 
fail to span X. 

The following experiment illustrates the behavior of FEAST for different numbers M. 
Experiment 3.1 

We consider a size-1059 matrix A = LAP_CIT_1059 from modelling cross-citations in scientific publications, 

and B = I. In this test we search for M = 1 , . . . , 450 eigenpairs with eigenvalues in an interval Ix containing 
the M = 300 lowest eigenvalues. The maximum number of iterations allowed for FEAST is 20. 

The left panel of Figure [T] shows the number of iterations necessary for FEAST to calculajte all eigenpairs 
within Ix with sufficiently small residual \\Ax — \Bx\\ < £-ri-max{|A| , |A|}, as a function of M. An iteration 
count of 20 typically implies that either none or not all eigenpairs converged within these 20 iterations. The 
right panel shows the residual span for all computed eigenpairs with eigenvalues in the interval after the 
respective number of iterations (20 or fewer, if convergence was reached beforehand). Again these numbers 
are given as a function of M . We see that, leaving aside the very small region arouncFthe exact eigenspace 
size, either all or none of the eigenpairs show a sufficiently sniall residual. While for M < M no eigenpairs 
converge and especially the minimum residuals are large, for M > M also the maximum residuals begin to 
drop significantly and typically all eigenpairs may converge if only enough iterations are performed. With 
M just slightly larger than Af , all eigenpairs reach convergence within few iterations. 

For a better understanding of the evolution of the computed eigenspace, we monitored the largest canoni- 
cal angle < (X^*^ X/^) p. 603] between the current approximate eigenspace X'^^^ and the exact eigenspace 
Xi^, as well as the angle < (X'^*^ X^*"-*^^) between the current and the previous iterate. Figure |2] provides 
these angles for three values of Af, M ~ 250, M — 300 and A/ — 350. In this last case, after five iterations 
the computed eigenspace contains tiie exact one and does not vary anymore; these two facts imply conver- 
gence. By contrast, the curves for M — 250 indicate that while the computed eigenspace becomes contained 
in the exact one after more than 20 iterations, it keeps varying, never to reach^ convergence. Interestingly, the 
worst convergence with respect to the exact eigenspace seems to occur for M = 300. This can be intuitively 
understood by the fact that two subspaces of the same dimension need to be identical in order to have an 
angle of zero between each other. 

3.2. Choice of the starting basis Y 

The choice of the starting basis Y e C"xAf j-^ line[l]of Algorithm [l] plays a critical role in the unfolding of the 
algorithm: most importantly, Y has to include components along U , so that its projection through Q spans 
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Figure 1: Left: Number of necessary iterations. Right: Minimal (lower line) and maximal (upper line) residual. 
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Figure 2: Canonical angles (left: between current iterate X'*' and exact eigenspace Xj; right: between current iterate X(') 
and previous iterate X^'"^)) in degrees for A/ = 250,300,350. 



X. Vice versa, if one or more of the columns of Y are i?-orthogonal to the space U, then the corresponding 
columns of [/ = QY will be zero. If a good initial guess for the eigenvectors of {A, B) is available, then it 
can used as the starting base Y] otherwise the typical choice is a base of random vectors. 

Experiment 3.2 

We used FEAST with M = 450 to compute the M = 300 lowest eigenvalues and corresponding eigenvectors 
of the matrix pencil {A = LAP_CIT_1059, S — I). These eigenvalues are simple and sufhciently far away 
from the only multiple eigenvalue, which is zero. With a fixed random starting basis Y , four iterations were 
sufficient to compute all wanted eigenpairs with residuals \\Axj — XjBxj\\ < 5.5 x 10""'^^. Then we projected 
out the ten eigenvectors corresponding to the 10 lowest eigenvalues via Y :— {I — Xi-ioX^.j^q) ■ Y. It took 
seven iterations for the lowest 290 eigenpairs to converge, and nine more iterations for other nine. One 
eigenpair did not converge within the limit of 20 iterations. Convergence took place "from top to bottom", 
i. e., the eigenpairs 11, ... , 300 converged first, then the eigenpairs 2, . . . , 10. The smallest eigenvalue did not 
converge within the iteration limit. 

In general, thanks to round-off errors, convergence could still be reached in most of our tests. In fact 
even though some components were zeroed out, the floating point arithmetic causes almost zero entries to 
grow as the computation progresses. Convergence is then reached with noticeably more iterations. 
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3.3. Stopping criteria 



Algorithm 1 relies on a stopping criterion to determine whether the eigenpairs are computed to a sufficient 
degree of accuracy. Such a criterion must balance cost and effectiveness, fn the original implementation PQ, 
FEAST monitors convergence through the change in the sum of the computed eigenvalues; more precisely, 
a relative criterion of the form 

Itracei, — tracei._i I , , 

' J < TOL (5) 

Itracefcl 

is used, where trace^ denotes the sum of the Ritz values in the fcth iteration lying in the search interval and 
TOL is a user tolerance. Criterion ^ raises three problems. 

First, its denominator might be zero or close to, causing severe numerical instabilities. Such a scenario 
arises, for instance, when all the eigenvalues in the search interval are zero. Second, the number tracer can 
be (almost) constant for two consecutive values of k, stopping FEAST even if the residuals are still large. 
The third problem is of a more general nature: if the algorithm stagnates before the considered eigenpairs 
converge, all criteria that are based only on the change in the eigeni)a/Mes might still signal convergence. It 
is not hard to construct examples where this happens 1^. 

As an alternative criterion we propose a per-eigenpair residual that depends on the search interval: a 
Ritz pair has converged if it fulfills the inequality 

llAx - SxAII < e • n -maxllAI , |A|}. (6) 

The extra cost of this criterion is one matrix-vector product per^ector because Bx is needed anyway to 
compute Y in the next iteration, or, if sparsity is not exploited, 0{M'^n) operations since Ax — {AU)w, with 
AU reused from Step 2 of the FEAST algorithm. If max{|A| , |A|} is too small, e.g., when only eigenpairs 
with zero eigenvalues are sought after, one may replace this quantity by an estimate for ||_B~^A||, which is 
the magnitude of the largest eigenvalue of the problem and might be obtained by some auxiliary routine 
(e.g., by some steps of a Lanczos method, see ^15]). 

Let us consider a matrix with spectrum symmetric with respect to zero, and choose the search interval 
to be symmetric around zero. If summed, the computed Ritz values cancel themselves pairwise, thus traccfe 
is approximately zero. When FEAST is applied to such a matrix, the trace criterion signals convergence 
when the difference in ^ comes close enough to zero; in this case, the ratio ^ holds no information about 
convergence. From multiple tests with random starting bases, we experienced that it took between 7 and 
100 iterations for criterion ([5| to signal convergence. By contrast, criterion ([6| dropped below 10^^^ always 
after 6 iterations. 

In an additional test, we ran FEAST on a symmetric matrix of size 470, ^eking the known 57-fold 
eigenvalue 1. Here we chose the search interval symmetrically around 1, and M = 120 > M. After five 
iterations, ([5| was 1.2 x 10~^^, effectively halting the computation, although the residuals (|6| were still of 
order 10~® to 10~^^. The right hand side of ^ was about 10"^'^ in this example, meaning that none of the 
residuals was satisfactory small. 

3.4. The impact of the hnear solver on residuals and orthogonality 

As seen in Section |2.3| the computation of the basis according to ([T]) involves solving several linear systems 
of the form 

{zB - A)V = BY (7) 

for V. The particular values of the integration points z depend on the method chosen for numerical inte- 
gration, which is not discussed here. Typically, z will be a complex number near the spectrum of {A,B). 
Recall that ([t]) is a linear system with AI right hand sides. In principle, any linear solver can be used. Direct 
solvers, e. g., Gaussian elimination based, can be prohibitively expensive because for each value of z we need 
to factorize zB — A, which may be an 0{n^) process. 

The methods of choice for solving large sparse linear systems without further knowledge about the 
underlying problem are Krylov subspace methods; for a review see, e.g., |10j . Here we cannot give a 
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detailed discussion, but let us remark that the convergence of Krylov subspace methods depends on several 
parameters. First, the best convergence results can be expected for Hermitian matrices since this property 
can be exploited. Unfortunately, ([t]) typically has a non-real diagonal and therefore is not a Hermitian 
problem. (However, if B = / then the matrix zl — A is shifted Hermitian, so methods for shifted systems 
may be applicable |14j.l Second, the convergence for a fixed method typically depends on the structure of 
the spectrum of the matrix. The eigenvalues of zB — A are scattered over the complex plane so that no good 
convergence results can be inferred. Third, the condition number \\izB — ^)^^|| • \\zB — A\\ of the system 
plays an important role and is often large for ^ , since z can be very close to the spectrum of {A, B) . For 
these reasons, standard Krylov subspace solvers may need a large number of iterations to converge. This 
expectation was confirmed by our experiments. The need for an effective preconditioner is apparent, and 
its development is part of further research. 

Another way to speed up the linear solvers is to terminate them before full convergence is reached. Thus 
the question arises how accurately the systems ([t]) need to be solved in order to obtain eigenpairs of sufficient 
quality in a reasonable number of FEAST iterations. We therefore investigated the effect of the accuracy in 
the solution of linear systems on the ultimately achievable per-eigenpair residuals and the orthogonality of 
the eigenvectors, as well as on the number of FEAST iterations. 

Experiment 3.3 

We applied Algorithm 1 to the matrix pair {A,B), where A = LAP_CIT_395 arises in the modelling of 
cross-citations in scientinc publications, and B was chosen to be a diagonal matrix with random entries. 
We calculated the eigenpairs corresponding to the 10 largest eigenvalues. The linear systems were solved 
column-by-column by running GMRES [TD] until \\{zB - A)vj - ByjW / ||r°|| < £ii„, where r° is the starting 
residual. Figure |3] reveals that the residual bounds that were required in the solution of the inner linear 
systems translated almost one-to-one into the residuals of the Ritz pairs. Even for a rather large bound 
such as e\in = 10^^, the FEAST algorithm still converged (even though to a quite large residual). For 
the orthogonality of the computed eigenvectors xj, the situation was different. After 20 FEAST iterations, 
an orthogonality level max^^j |a:Pi3xj| of order 10""'^^ could be reached for each of the bounds Sun = 
10~^, 10"®, 10~^°, 10"^^ in the solution of the linear systems. Thus the achievable orthogonality does not 
seem to be very sensitive to the accuracy of the linear solves. It also did not deteriorate significantly for a 
larger number of desired eigenpairs. 




Figure 3: Range of all residuals among all Ritz pairs in Ix for four different residual bounds eu^ in the linear solver. 
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4. FEAST for multiple intervals 



Using FEAST for computing a large number M of eigenpairs is not recommended since the complexity grows 
at least as 0{nA'P) due to the matrix-matrix products in Steps 2 and 4 of Algorithm[l] or 0{n^M) if sparsity 
is not exploited, because 0{M) matrix-vector products must be computed. (In fact, if M approaches n 
then just Step 3 of Algorithm [T] has roughly the same complexity as the computation of the full eigensystem 
of the original problem.) However, as already hinted at in |T], FEAST's ability to determine the eigenpairs 
in a specified interval makes it an attractive building block to compute the eigenpairs by subdividing the 
search interval I\ into K subintervals and applying the algorithm — possibly in parallel — to 

each one of them. 

In the following we report on issues concerning the orthogonality of eigenvectors coming from different 
subintervals. It is well known that eigenvectors computed independently from each other tend to have worse 
orthogonality than those obtained in a block-wise manner. Furthermore, the quality of the results depends 
on the internal structure of the spectrum, namely the relative distances between the eigenvalues. For more 
details, see, e.g., [T5] . 

In the following we distinguish between global and local orthogonality, according to the definitions 
orthgiobai = max \xf BxA and orth/; = max |a;f^-Bxj| . 

Note that orthgiobai denotes the worst orthogonality among all computed eigenvectors, while orth^ describes 

(k) 

the orthogonality achieved locally for the fcth subinterval I\ . The next two experiments reveal quite 
different behavior, depending on the presence of clusters and on the choice of subintervals. 

Experiment 4.1 

In this experiment we calculate the 800 lowest eigenpairs of the size-1473 matrix pair (bcsstkll, bcsstmll) 
from the Matrix Market (http://math.nist.gov/MatrixMarket/). The corresponding eigenvalues range 
from 10.5 to 3.8 x 10^ and are not clustered. We utilize different numbers of subintervals, K = 1, . . . , 5, 
and K = IQ. Figure |4] shows that while the local orthogonality is high and maintained as the number of 
intervals increases, the global one degrades by two or more orders of magnitude. 




Number of subintervals 

Figure 4: Global orthogonality and range of local orthogonality orthj for K = 1, . . . ,b and K = 10. 

Experiment 4.2 

In this test we consider a real unreduced tridiagonal matrix A of size 2003. Its eigenvalues are simple, 
even though some are tightly clustered; see the top plots in Figure |5] The objective is to compute the 300 
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Figure 5: Computation of the eigenpairs corresponding to the 300 largest eigenvalues A1704, . . . , A2003 with the subdivision 
point /i = Aig25 taken from a group of very close eigenvalues. The top plots show the eigenvalues and the subdivision points 
(vertical lines), the bottom plots give a pictorial visualization of the orthogonality x^Bxj , i j^l. 



largest eigenpairs. To this end we initially split the interval I\ = [A1704, A2003] into /| = [Ai704,/i] and 
I^^-' = A2003]) with /X = Ai825 ~ 0.448 x 10~^ chosen within a a cluster of 99 eigenvalues. The relative 
gap between eigenvalue A1825 and its neighbors is of about 10^^^ (i-e., agreement to roughly eleven leading 
decimal digits) . A sketch of the eigenspectrum with /i is given in the top left of Figure [s] 

While FEAST attains very good local orthogonality for both subintervals (orthi — 4.4 x 10^^^ and 
orth2 ~ 5.7 X 10"^'*), it fails to deliver global orthogonality (4.7 x 10^*). In the bottom left plot of Figure [s] 
we provide a pictorial description of |a::^_Ba;j|, Xi,Xj G Ix- The dark colored regions indicate that the loss 
of orthogonality emerges exclusively from eigenvectors belonging to the 99-fold cluster. Next we divide the 
interval into 3 segments making sure not to break existing clusters (see top right of Figure |5|. As illustrated 
in the bottom right plot, both the local and global orthogonality are satisfactory (10~^^ or better). 



5. Conclusions 

We expounded the close connection between FEAST algorithm and the well-established Rayleigh-Ritz 
method for computing selected eigenpairs of a generalized eigenproblem. Starting from the mathematical 
foundation of this connection, we identified aspects of the solver that might play a critical role for its accuracy 
and reliability. Specifically, we discussed the choice of starting basis and stopping criterion, and the relation 
between the accuracy of the solutions of the linear systems internal to FEAST and the resulting eigenpairs; 
we also investigated the use of FEAST for computing a large portion or even the entire eigenspectrum. 
Through numerical examples we illustrated how each of these aspects might affect the robustness of the 
algorithm or diminish the quality of the computed eigensystem. 

While we hinted at possible improvements for several of the existing issues, some questions remain open 
and are the subject of further research. For instance, a mechanism is needed to overcome the problem 
of having to specify both the boundaries of the search interval and the number of eigenvalues expected. 
Additionally, it would be desirable to have a flag assessing whether all the eigenpairs in the search interval 
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have been found. Finally, orthogonality across multiple intervals should be guaranteed. 

In summary, our findings suggest that at the moment FEAST is a promising eigensolver for a certain 
class of problems, i.e., when a small portion of the spectrum is sought and knowledge of the eigenvalue 
distribution is available. On the other hand, we believe it is still not yet competitive as a robust "black box", 
general-purpose solvei[^ 
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^At the time of submission, v2.0 of the FEAST Package Solver was announced. This new version is designed for parallel 
computation, and doesn't seem to address the issues raised in this research paper. We still would like to remark that all the 
numerical experiments in this paper were executed making use of the algorithm from FEAST vl.O. 
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